!DIR$ ATTRIBUTES FORCEINLINE :: pfqsdt
  pure real(rkx) function pfqsdt(t,p) result(qsdt)
!$acc routine seq
    implicit none
    real(rkx) , intent(in) :: t             ! Temperature (K)
    real(rkx) , intent(in) :: p             ! Pressure (Pa)
    real(rkx) :: td, es, esdt, vp, vp1, vp2
    !
    ! For liquid water (temperature range 0C-100C)
    !
    real(rk8), parameter :: a0 =  0.611213476e+03_rk8
    real(rk8), parameter :: a1 =  0.444007856e+02_rk8
    real(rk8), parameter :: a2 =  0.143064234e+01_rk8
    real(rk8), parameter :: a3 =  0.264461437e-01_rk8
    real(rk8), parameter :: a4 =  0.305903558e-03_rk8
    real(rk8), parameter :: a5 =  0.196237241e-05_rk8
    real(rk8), parameter :: a6 =  0.892344772e-08_rk8
    real(rk8), parameter :: a7 = -0.373208410e-10_rk8
    real(rk8), parameter :: a8 =  0.209339997e-13_rk8
    real(rk8), parameter :: b0 =  0.444017302e+02_rk8
    real(rk8), parameter :: b1 =  0.286064092e+01_rk8
    real(rk8), parameter :: b2 =  0.794683137e-01_rk8
    real(rk8), parameter :: b3 =  0.121211669e-02_rk8
    real(rk8), parameter :: b4 =  0.103354611e-04_rk8
    real(rk8), parameter :: b5 =  0.404125005e-07_rk8
    real(rk8), parameter :: b6 = -0.788037859e-10_rk8
    real(rk8), parameter :: b7 = -0.114596802e-11_rk8
    real(rk8), parameter :: b8 =  0.381294516e-14_rk8
    !
    ! For ice (temperature range -75C-0C)
    !
    real(rk8), parameter :: c0 =  0.611123516e+03_rk8
    real(rk8), parameter :: c1 =  0.503109514e+02_rk8
    real(rk8), parameter :: c2 =  0.188369801e+01_rk8
    real(rk8), parameter :: c3 =  0.420547422e-01_rk8
    real(rk8), parameter :: c4 =  0.614396778e-03_rk8
    real(rk8), parameter :: c5 =  0.602780717e-05_rk8
    real(rk8), parameter :: c6 =  0.387940929e-07_rk8
    real(rk8), parameter :: c7 =  0.149436277e-09_rk8
    real(rk8), parameter :: c8 =  0.262655803e-12_rk8
    real(rk8), parameter :: d0 =  0.503277922e+02_rk8
    real(rk8), parameter :: d1 =  0.377289173e+01_rk8
    real(rk8), parameter :: d2 =  0.126801703e+00_rk8
    real(rk8), parameter :: d3 =  0.249468427e-02_rk8
    real(rk8), parameter :: d4 =  0.313703411e-04_rk8
    real(rk8), parameter :: d5 =  0.257180651e-06_rk8
    real(rk8), parameter :: d6 =  0.133268878e-08_rk8
    real(rk8), parameter :: d7 =  0.394116744e-11_rk8
    real(rk8), parameter :: d8 =  0.498070196e-14_rk8

    td = min(max(t - tzero,-75.0_rk8),100.0_rk8)
    if ( td >= 0.0_rk8 ) then
      es = min(a0 + td*(a1 + td*(a2 + td*(a3 + td*(a4 &
         + td*(a5 + td*(a6 + td*(a7 + td*a8))))))),0.15_rk8*p)
      esdt = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4 &
           + td*(b5 + td*(b6 + td*(b7 + td*b8)))))))
    else
      es = min(c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 &
         + td*(c5 + td*(c6 + td*(c7 + td*c8))))))),0.15_rk8*p)
      esdt = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4 &
           + td*(d5 + td*(d6 + td*(d7 + td*d8)))))))
    end if

    vp  = 1.0_rkx / (p - 0.378_rkx*es)
    vp1 = ep2 * vp
    vp2 = vp1 * vp
    qsdt = esdt * vp2 * p ! 1 / K
  end function pfqsdt

! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
